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CP ' We study the importance of hydrodynamic effects on the evolution of coalescing 

binary neutron stars. Using an approximate energy functional constructed 



o 



. from equilibrium solutions for polytropic binary configurations, we incorporate 
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hydrodynamic effects into the calculation of the orbital decay driven by gravitational 
wave emission. In particular, we follow the transition between the quasi-static, 
secular decay of the orbit at large separation and the rapid dynamical evolution of 
}^ • configurations approaching contact. We show that a purely Newtonian hydrodynamic 

c/3 . instability can significantly accelerate the coalescence at small separation. Such an 

instability occurs in all close binary configurations containing sufficiently incompressible 



stars. Calculations are performed for various neutron star masses, radii, and spins. 
The influence of the stiffness of the equation of state is also explored by varying the 
I effective polytropic index. Typically, we find that the radial infall velocity just prior 

to contact is about 10% of the tangential orbital velocity. Once the stability limit 
is reached, the final evolution only takes another orbit. Post-Newtonian effects can 
move the stability limit to a larger binary separation, and may induce an even larger 
radial velocity. We also consider the possibility of mass transfer from one neutron 
star to the other. We show that stable mass transfer is impossible except when the 
mass of one of the components is very small (M ^ 0.4Mq) and the viscosity is high 
enough to maintain corotation. Otherwise, either the two stars come into contact or 
the dynamical instability sets in before a Roche limit can be reached. 
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1. INTRODUCTION 



Coalescing neutron star binaries have long been recognized as a most promising source of 
gravitational radiation that could be detected by the new generation of laser interferometers 
such as LIGO (Clark 1979; Thorne 1987; Abramovici et al. 1992; Cutler et al. 1992). Statistical 
arguments based on the observed local population of binary pulsars and type lb supernovae 
lead to an estimate of the rate of neutron star binary coalescence in the Universe of about 
10~'^yr-i Mpc~3 (Narayan, Piran & Shemi 1991; Phinney 1991). Finn & Chernoff (1993) estimate 
that an advanced LIGO detector could observe about 70 events per year. In addition to providing 
a major confirmation of Einstein's theory of general relativity, the detection of gravitational waves 
from coalescing binaries at cosmological distances could provide the first accurate measurement of 
the Universe's Hubble constant and mean density (Schutz 1986; Cutler et al. 1992). Coalescing 
neutron stars may also produce intense bursts of neutrinos and high energy photons, and could be 
the source of extra-galactic gamma-ray bursts (Eichler et al. 1988; Paczyhski 1986, 1991). 

Recent calculations of the gravitational radiation waveforms from coalescing neutron star 
binaries have focused on the signal emitted during the last few thousand orbits, as the frequency 
sweeps upward from about 10 Hz to 1000 Hz. The waveforms in this regime can be calculated 
fairly accurately by performing high-order post-Newtonian expansions of the equations of motion 
for two point masses (Lincoln & Will 1990; Junker & Schafer 1992). Accuracy is essential here 
since the observed signal will be matched against theoretical templates. Since the templates must 
cover > 10^ orbits, a fractional error as small as 10~^ can prevent detection. When, at the end 
of the inspiral, the binary separation becomes comparable to the stellar radii, hydrodynamic 
effects become important and the character of the waveforms will be different. Special purpose 
narrow-band detectors that can sweep up frequency in real time will be used to try to catch the 
final few cycles of gravitational radiation (Meers 1988; Strain &; Meers 1991). In this terminal 
phase of the coalescence, the waveforms should contain information not just about the effects 
of relativity, but also about the internal structure of neutron stars. Since the masses and spins 
of the two stars, as well as the orbital parameters can be determined very accurately from the 
lower-frequency inspiral waveform, a simple determination of the maximum frequency fmax 
reached by the signal should provide a measurement of the neutron star radii and place severe 
constraints on nuclear equations of state (Cutler et al. 1993). Such a measurement requires 
theoretical knowledge about all relevant hydrodynamic processes. 
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Two regimes can be distinguished in which different hydrodynamic processes take place. The 
first regime corresponds to the 10 or so orbits preceding the moment when the surfaces of the two 
stars first come into contact. In this regime, the two stars are stih approaching each other in a 
quasi-static manner, but the tidal effects are very large. The second regime corresponds to the 
subsequent merging of the two stars into a single object. This involves very large departures from 
hydrostatic equilibrium, including mass shedding and shocks, and can be studied only by means 
of three-dimensional, fully hydrodynamic computations (Nakamura & Oohara 1991 and references 
therein; Rasio & Shapiro 1992) By contrast, in the first regime, the evolution of the system can 
be described fairly accurately by a sequence of near-equilibrium fluid configurations. We adopt 
such a description in this paper. Since neutron stars are not very compressible, the equilibrium 
conflgurations are not very centrally condensed and the usual Roche model for close binaries (e.g., 
Kopal 1959) does not apply. Instead, the hydrostatic equilibrium equation must be solved in three 
dimensions for the structure of the system. 

In two recent papers (Lai, Rasio & Shapiro 1993b, c, hereafter LRSl and LSR2), we have 
used an approximate energy variational method to study analytically the hydrostatic equilibrium 
and stability properties of Newtonian binary systems obeying polytropic equations of state. In 
particular, we have constructed the compressible generalizations of all the classical solutions for 
binaries containing incompressible ellipsoidal components (Chandrasekhar 1969). The so-called 
Darwin- Riemann configurations are of special importance for modeling neutron star binaries. 
They are generalizations of the Roche-Riemann configurations originally introduced by Aizenman 
(1968) to describe an incompressible star in circular orbit around a point mas^. Kochanek 
(1992) first pointed out the importance of these configurations for modeling neutron star binaries. 
In Darwin-Riemann binaries of the type considered here, fluid motions with uniform vorticity 
parallel to the rotation axis are present in the corotating frame of the binary system. These 
configurations represent a simple model of a close binary containing stars whose spins are not 
necessarily synchronized with the orbital motion. Indeed, it has been argued recently that the 
synchronization time of neutron star binaries may remain very long compared to their orbital decay 
timescale (Bildsten & Cutler 1992; Kochanek 1992). In the limit where the viscosity is negligible 
(leading to an infinite synchronization time), the fluid circulation is strictly conserved and an 
asynchronous conflguration of the Darwin-Riemann type is expected (Miller 1974; Kochanek 
1992). The value of the circulation is determined by the initial spins of the two neutron stars at 
large binary separation. If the stars spin uniformly at large separation, the vorticity will remain 
uniform. Similarly, if the spins are initially aligned with the rotation axis of the binary (a most 
likely configuration), the vorticity vector will remain aligned as well. Detailed calculations of the 
properties of Darwin-Riemann binaries are presented in LRSl (for two identical components) and 
LRS2 (for nonidentical components). 



Roche-Riemann binary consists of a star and point mass, while a Darwin-Riemann binary consists of two 
finite-size stars. 
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An important result found in LRSl is that Darwin-Riemann configurations containing a 
sufficiently incompressible fluid can become hydrodynamdcally unstable. Close binaries containing 
neutron stars with stiff equations of state (T ^ 2) should therefore be particularly susceptible 
to these instabilities. As the dynamical stability limit is approached, the secular orbital decay 
driven by gravitational wave emission can be dramatically accelerated. It is this approach towards 
instability that we intend to study here. When the stability limit is reached, the two stars plunge 
almost radially toward each other, and merge together into a single object after just a few orbits. 
The complicated three-dimensional hydrodynamics characterizing this very short terminal phase 
of the evolution can only be studied with large-scale computer simulations (Nakamura & Oohara 
1991; Rasio & Shapiro 1992). In contrast, the long transition between the very slow secular decay 
of the orbit at large separation and the more rapid inspiral as the stability limit is approached is 
better explored using quasi-analytic methods. The analytic results should also prove helpful in the 
construction of better initial conditions for the numerical simulations. Indeed, nonsynchronized 
initial conditions are particularly difficult to deal with numerically (for recent attempts, see 
Shibata et al. 1993; Davies et al. 1993). 

In a previous study (Lai, Rasio, & Shapiro 1993a, hereafter Paper I), we modeled the 
hydrodynamic and tidal effects associated with finite-size, spinning components by modifying 
the orbital evolution equations for two point masses in a simple, heuristic way. In this paper, 
we improve our treatment by introducing the full binary equilibrium solutions for polytropes 
calculated in LRSl. Specifically, we model a neutron star binary as a Darwin-Riemann 
configuration containing two triaxial polytropes. We calculate the orbital evolution of the system 
as the stability limit is approached and we study how the radial infall evolves from secular 
to dynamical. Because of the simplicity of our method, we can perform a large number of 
calculations, systematically surveying all relevant parameters such as the stellar masses, radii, and 
spins. We also vary the polytropic index of the fluid to mimic the mean stiffness of the nuclear 
equation of state. All of these parameters can be chosen separately for each of the two companion 
stars. However, we focus our attention on binaries containing two identical neutron stars. Indeed, 
all the measured masses of neutron stars in our Galaxy appear to be narrowly clustered around 
1.4M0 (Thorsctt et al. 1993). An important simplifying assumption which remains in the present 
treatment is that at any separation the internal degrees of freedom of the two stars (surface 
deformations and central densities) are calculated from the values in circular equilibrium while we 
follow the orbital decay. This approximation is reasonable whenever the orbital timescale remains 
shorter than the radial infall timescale, but longer than the internal dynamical timescale of each 
star. In a future paper, we will derive and solve the generalized evolution equations for a binary 
configuration in which all degrees of freedom are allowed to vary dynamically. 

In §2 we present our Darwin-Riemann equilibrium model for binaries containing two identical 
stars. We also mention briefly Darwin-Riemann models for two stars with unequal masses, as 
well as Roche-Riemann models of neutron-star-black-hole binaries. The orbital evolution of these 
models is calculated in §3 and §4. General relativistic corrections to the orbital motion are 
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discussed in §5. In §6, we examine the different types of terminal evolutions that are possible 

when the two stars have different masses. In particular, we address the possibility that stable 
mass transfer from one neutron star to the other may occur (as first suggested by Clark &; Eardley 
1977). 
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2. DARWIN-RIEMANN EQUILIBRIUM MODELS FOR CLOSE BINARIES 



In this section, we derive equilibrium equations for compressible Darwin-Riemann binary 
configurations containing two identical polytropes. We use an energy variational method together 
with an ellipsoidal approximation to obtain those equations. Under the combined effects of 
centrifugal and tidal forces, each polytrope assumes a nonsphcrical equilibrium configuration 
which we approximate by a triaxial ellipsoid. As in the classical Darwin problem, we assume 
that the two ellipsoidal figures corotate with the orbital angular velocity (i.e., we consider only 
stationary tides), but, as in the classical Riemann ellipsoids of type S, we allow for internal fluid 
motions with uniform vorticity parallel to the rotation axis. More details about the method and a 
number of applications to other types of binary and isolated rotating configurations can be found 
in LRSl and LRS2. 

Consider a binary system containing two identical polytropes of mass M and polytropic index 
n, in circular orbit about each other. The density and pressure are related by 



where the adiabatic exponent F = 1 + 1/n. The constant K measures the specific entropy and is 
the same for both stars. We denote the central density by pc and the binary separation by r. We 
make the assumption that the surfaces of constant density within each star can be approximated 
by self-similar ellipsoids. The geometry of the configuration is then completely specified by the 
three principal axes of the outer surfaces, ai, 02, and 03, with ai measured along the axis of the 
binary, 02 in the direction of the orbital motion, and 03 along the rotation axis. In addition, we 
assume that the density profile p{m), where m is the mass interior to an isodensity surface, is 
identical to that of a spherical polytrope of same volume and entropy. 

Under these assumptions, the total internal energy in the system is given by 



2.1. Basic Equations 



P = Kp^, 



(1) 




(2) 



and the total self-gravitational potential energy (setting G = 1) is 




(3) 




(4) 
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such that / = 1 for spherical stars. The dimensionless coefficients Ai are defined as in 
Chandrasekhar (1969, §17); they depend only on the axis ratios a^/ai and 03/02- The coefficients 
ki and k2 are dimensionless polytropic structure constants which depend only on n, 

n(n+l) 3 f 4n\e\\ \'/' 

where 9 and are the usual Lane-Emden variables for a polytrope (see, e.g., Chandrasekhar 1939). 
The gravitational interaction energy Wi between the two stars is given, to quadrupole order, 

by 

M 

VF, = -— -^(2/11-/22-/33), (6) 

where 

hj = Kn—^Oij (no sum) (7) 

and constant depending only on n, 

'-' ' -in M 



3^1 |y i| Jo 

so that K„ = 1 for n = 0. Values of ki, ^2 and k„ for different n are given in Table 1 of LRSl. 

Now turn to the kinetic energy. For a synchronized binary, the fluid is simply in uniform 
rotation at the orbital frequency i}e^ perpendicular to the orbital plane. Here, however, we allow 
for an additional internal motion of the fluid inside each star, assumed to have uniform vorticity 
Ce3 as measured in the corotating frame of the binary, 

C^(Vxu).e3 = -^i±^A, (9) 
0102 

where 

01 02 . 

u = — Ax2ei Axie2 (10) 

02 ai 

is the fluid velocity in the corotating frame. Here the origin of the coordinates Xj is at the center 
of mass of the star. The quantity A is the angular frequency of the internal fluid motions. Note 
that the velocity fleld ( p!o| ) is everywhere tangent to the ellipsoidal surfaces of constant density. 
For a synchronized binary system, we have u = = A = 0. The velocity fleld in the inertial frame 
relative to the center of mass of the star is given by 

u(°) = u + nxx. (11) 

The vorticity in the inertial frame is 

C(°)^(VxuW).e3 = (2 + /R)0, (12) 

where we have deflned the ratio 

fn - ^. (13) 
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It is straightforward to calculate the total kinetic energy corresponding to this velocity field. We 
find 



T = Ts + To = 2 



+ i^r^l^^, (14) 



2' 

where fj, = M/2 is the reduced mass and / = In + I22 = KnM{ai + a^/S is the moment of inertia 
of each star. We have also defined the orbital kinetic energy To = (^/2)r2j72 and the "spin" kinetic 
energy = T — Tq. The total angular momentum J can be written similarly as 

2 

J = 2{m K„Maia2A) + /ir^O. (15) 

5 

Another important conserved quantity is the fluid circulation C along the equators of the two 
star. Following LRSl we write 

C = (-^KnM\c = 2f--^K„MVaia2C^°^ = 2(/A - ^K„Maia2l^). (16) 
V ovr / V OTT / 5 

Note that C has dimensions of an angular momentum but is proportional to the conserved 
circulation. We usually refer to C itself as the circulation. Using equations (|l5|) and (|l^), the 
kinetic energy can be rewritten as 

T = T, + To = T+ + r_ + ^^r2f]^ (17) 

where 

T± = i/±(l^ ± Kf = ^( J ± C - ixr'^nf, (18) 



and 



/± = l^KnMiai ^ a2f = 2Is/h±, (19) 
5 



with Is = 2K„Mi?2/5 and h± = 2R^/{ai T 02)^- 

The total energy of the system (not necessarily in equilibrium) can now be written simply as 
the sum of expressions @, (H) , (P) , and dl] 



^(Pc,Ai,A2,r; M, J,C) = U + W + Wi + T, (20) 

where we have chosen as independent variables the central density p^, the binary separation r, and 
the two oblateness parameters Ai = (03/01)^/^ and A2 = (a3/a2)^''^- The equilibrium structure of 
the binary can be determined from the four conditions 

dr dpc dXi 8X2 

The first condition, dE/dr = 0, gives the equilibrium relation between and r, i.e., the modified 
Kepler's law for the binary: 

o 2M 

= —{1 + 25) =2pr{1 + 25), (22) 
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where we have defined 



3 (2/11 - I22 - hs) _ ,1^/ 3 



The second condition, dE/dpc = 0, leads to the virial relation for the binary, 



-U + W + 2Ts 
n 



2M^ 



R 



-9t, 



where we have defined 



R 



9t = 



2 R 



11 — ^22 — -133 J 



5. 



From equations (^), (|3|) and (p^), the equilibrium mean radius R can be obtained as 

-n/(3— n) 



R — Rq 



1-2^1/ 



5 — n 



gt 



\W\J' V 3 

where Rq is the radius of a spherical equilibrium polytrope with the same mass and entropy. 



(23) 

(24) 
(25) 

(26) 



-(l-n)/(3-n) 



{n+l)K 
4tt 



n/{3~n) (l-n)/(3-n) 

47r, 



(27) 



The last two conditions dE/d\i = dE/d\2 = 0, together with the virial relation, can be used 
to derive two equations determining the axes ratios in equilibrium. 



where we have defined 



—al + 2 ( 2 + 25 + 



Q2^l \ 2 , 2 
ai + 03 



^a?+(l + 4<5 



Qi 



2Q1Q, \ 2 I 2 

^2 + «3 

fJ-R 



2{alAi-alA^), 
2ialA2-alA3), 



al 



ai 



2 , 2^ = +— A, 
af + 02 a2 



al 



Q2 = +^ , 



ai 



(28) 
(29) 

(30) 
(31) 



and Qn = Kn{l - n/5), fiR = ^ir/{'kp), with p = 3M/(47raia2a3). 

For any given M, J, C, K, and n, an equilibrium model can be constructed from the four 
algebraic equations (p2[), (26), ( |28| ) and (p9|). For specified values of /r and r/oi, we solve 
equations ( [2^ ) and ( [2^ ) for the axis ratios 02/01 and 03/01 by a Newton- Raphson method following 
an initial guess. The mean radius of the star is then obtained from equation (|2^). The total 
equilibrium energy of the system can be written 



1 



„2r.2 



Eeq = 2Es + -pr^n 



M2 
r 



/ 2n + 3\ M 



-{21 



11 — -^22 



'33J, 



(32) 
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where Eg is the intrinsic energy of each star, 

(3-n)M2 



f 



2n 



n 



W 



{5-n)R 

The total equihbrium angular momentum J^g is given by equation (^), evaluated for the 
equilibrium solution. 



(33) 



2.2. Equilibrium Sequences with Constant C 



Since the circulation C is conserved in the absence of viscosity, it is useful to construct 
sequences of equilibrium configurations with fixed C. A constant-C sequence is parametrized by 
the binary separation r or the angular momentum J. Note that in general fn varies along a 
constant-C sequence. When C is specified, the value of fn needs to be determined simultaneously 
with 02/01 and 03/01 to satisfy equation ([T^). 



Of particular interest here is the irrotational Darwin-Riemann sequence, for which the 
circulation C = (//j = —2). This corresponds to the case where the stars have no spin at large 
separation. In Table 1 we have listed some properties of these irrotational Darwin-Riemann 
configurations for n = 0, 0.5, 1, and 1.5. All sequences are terminated when the surfaces of the two 
stars are in contact, i.e., r/oi = 2. Following LRSl, we adopt units based on the radius Rq of a 
spherical equilibrium polytrope with same mass and entropy (eq. [^), defining the dimensionless 
quantities 

- Q - J - E 

^^VpoF"' '^^WRo)^' ^^WJRo)' ^^^^ 

where po = M/{ATrRl/3). Note that R/Ro > 1 for n > 0, indicating that the volume of each star 
increases when placed in a binary system. 

Equilibrium sequences with constant C 7^ can also be constructed. The value of C can be 
specified from the spin angular frequency of each star at large separation. Indeed, for large r, 
we have 01^02,^^^^ 2M/r^ and 

J fir'^n-2IA = nr'^n + 2ms, (35) 
C 21 A = -2ms, (36) 

where we have identified fi^ = — A(r = 00) as the spin angular velocity at large r (for an 
axisymmetric star, uniform spin and vorticity are indistinguishable in the ellipsoidal models). Note 
that when fi^ is positive (i.e., the spin is in the same direction as the orbital angular momentum), 
C is negative. The maximum spin that a uniformly rotating neutron star can sustain without 
shedding mass from its equator is given by (Friedman, Ipser & Parker 1986; Cook, Shapiro & 
Teukolsky 1992) 
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For typical neutron stars, this maximum spin rate is not very sensitive to the adopted equation 
of state. Table 2 gives the equilibrium properties of the constant-C Darwin-Riemann sequences 
corresponding to 0^ = 0.2 and 0.4, for a polytropic index n = 0.5 or 1. Irrotational cases, 
corresponding to f^^ = 0, were given in Table 1. 

In Figure 1, we show the variation of the total equilibrium energy Eeq{r) and angular 
momentum Jeqif) of the binary system, as well as the orbital angular frequency O, along several 
sequences with constant C = —2IQ,s, for Clg = 0, 0.1, 0.2, 0.4 and polytropic index n = 0.5. 
Figure 2 shows the variation of the principal axes along the same sequences. In these plots, each 
curve terminates at the contact solution. We see immediately that there exists a critical separation 
rm where E^qir) and Jeq{r) are simultaneously minimum. That Egq and Jeq attain their minimum 
simultaneously is a consequence of the property dEeq = 0, d Jgg along an equilibrium sequence with 
constant C (LRSl, Appendix D). The minimum occurs as a result of the strong tidal interaction 
between the two stars at small separation (see Paper 1 for a qualitative discussion). 

The minima in Eeq{r) and Jeq{r) along a constant-C sequence indicate the onset of dynamical 
instability. Indeed, at r = r^, it becomes possible for a small dynamical perturbation of the 
system (which conserves C) to cause no first-order change in the equilibrium energy or angular 
momentum. Such a perturbation must have eigenfrequency a;^ = 0, signaling the onset of 
instability (see, e.g., Shapiro &: Tcukolsky 1983, Chap. 6; Tassoul 1978). More rigorously, it can 
be shown (LRSl) that the onset of instability, determined from the condition 

(d'^E \ 
— — - — I =0, i, J = 1, 2, . . . (onset of instability), (38) 
daidaj ) 

where the aj's are the parameters specifying the configuration (in this case, r,Pc,Xi and A2), 
exactly coincides with the points of minimum Eeq{r) and Jeq{r). Binary configurations with 
r < rm are thus dynamically unstable. From Tabic 1, wc see that is smaller for larger n. This 
is because tidal effects become important at smaller separation for more centrally-concentrated 
stars. When n ^ 1.2, the minimum disappears and all binary configurations with C = remain 
stable up to contact. For sequences with C 7^ 0, we see from Table 2 that, for a given n, the 
minima become more shallow and ultimately disappear as \C\ increases. Qualitatively, this is 
because the intrinsic oblateness of the two spinning stars causes them to come into contact at a 
larger separation, where tidal effects are smaller. 



2.3. Synchronized Equilibrium Sequences 

When viscosity is important, the circulation C is no longer conserved during the evolution of 
the binary. In the limit where the synchronization timescale is much smaller than the orbital decay 
timescale, the evolution of the system may be described approximately by a sequence of uniformly 
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rotating (i.e., synchronized) equilibrium configurations. This is the case for the vast majority of 
observed binaries, but may not be true for neutron star binaries (see Bildsten & Cutler 1992; 
Kochanek 1992). Nevertheless, as a limiting case, we also construct equilibrium sequences with 
uniform rotation. These represent the compressible generalizations of the classical incompressible 
Darwin configurations discussed in Chandrasekhar (1969). The corotating configurations can be 
constructed as a special case of the Darwin-Riemann solutions with fn = and A = 0. The total 



kinetic energy (eq. [17|) in this case simply reduces to 



T=- (/ir^ + 2I)Q'^ = —4^ . (39) 

2^^ ' 2(/xr2 + 2I) ^ ' 

Extensive tables and plots illustrating the properties of compressible Darwin sequences are given 
in LRSl and will not be repeated here. Note that the circulation C varies along those sequences. 

As in Figure 1, it is found that E(,q{r) and Jeqif) can attain a minimum at a critical 
separation r^ before contact is reached. This occurs for sufficiently incompressible configurations, 
with n ^2. The minimum is more pronounced in this case because of large positive contributions 
to the Egg and J^q from the synchronized spins. 

Here the minimum marks the onset of secular instability along the sequence. In the presence 
of viscosity, configurations with r < r^ will be driven away from synchronization (see Counselman 
1973, and Hut 1980, for simple models of secularly unstable binaries). The instability at 
T = rm cannot be dynamical because neighboring configurations along the sequence are still in 
uniform rotation and therefore can only be reached on a viscous timescale (recall that dynamical 
perturbations conserve C, which varies along the corotating sequence). True dynamical instability 
occurs a little further along the sequence (at slightly smaller r), when neighboring configurations 
with the same value of C can be reached with no change in equilibrium energy to first order (see 
LRSl, in particular Fig. 14). In this paper (as in Paper 1), for simplicity, we do not distinguish 
between the secular and dynamical stability limits and we treat the instability at r = as if it 
were dynamical. This is justified because the binary separation changes very little between the 
secular and dynamical stability limits, and departures from synchronization remain always small. 
See LRS2, however, for a detailed discussion of what really happens between the secular and 
dynamical stability limits. 



2.4. Models of Black Hole — Neutron Star Binaries 



Binary systems consisting of a point-like object orbiting a finite-size star, such as a black 
hole - neutron star (BH-NS) binary, can be modeled as Roche-Riemann configurations. We have 
studied such configurations in detail in LRSl. Compared to the equal-mass Darwin-Riemann 
configurations discussed in §2.1-2.3, a prominent new feature in the Roche-Riemann systems is the 
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existence of a Roche limit prior to contact for circular equilibrium. The Roche limit corresponds 
to the point at which the binary separation r has a minimum value below which no equilibrium 
solution exists. At the Roche limit, the slope of the E(,q[r) and Jeq{r) curves becomes infinite. 
Beyond the Roche limit, there exists a second branch of equilibrium solutions, with larger surface 
deformation, that extends all the way to contact. However, these solutions beyond the Roche limit 
are unphysical since they have higher energies than those along the main equilibrium branch for 
the same value of r (cf. Figs. 10, 11, and 13 of LRSl). When the orbit of a binary decays to the 
Roche limit, tidal disruption or mass transfer is unavoidable. 

For Roche-Riemann equilibrium sequences with constant C, a dynamical instability is always 
encountered prior to the Roche limit. Indeed, the dynamical stability limit corresponds to the 
minimum in the E(.q{r) curve, which must always precede the point of infinite slope. Thus 
binaries at the Roche limit are always dynamically unstable. In contrast, along a sequence 
of synchronized configurations (a compressible Roche sequence), the minimum of Eeq{r) only 
corresponds to a secular instability, as discussed in §2.3. Dynamical instability sets in later (at 
smaller r) and the Roche limit may or may not be dynamically unstable, depending on the mass 
ratio and the compressibility of the star. For typical BH-NS binaries with moderate mass ratios 
{Mbh /Mns — 10) and highly incompressible neutron star matter (n ~ 0.5), we find that the 
dynamical instability sets in prior to the Roche limit (see Table 10 of LRSl). Consequently, these 
BH-NS binaries are always dynamically unstable at the Roche limit. Only when the neutron star 
is very compressible or when the black hole is much more massive, can the Roche limit set in prior 
to the dynamical stability limit. In those dynamically stable binary at the Roche limit 

can exist. 

In the rest of this paper, we focus on NS-NS binaries, rather than BH-NS binaries, for the 
following reason. The onset of instability (secular or dynamical) along a Roche or Roche-Riemann 
sequence occurs at ~ 2(1 + qy^^R, where q = Mbh /^NS- But the last stable circular orbit 
around a Schwarzschild black hole is at rcR ~ 6Mbh for Mbh ^ ^NS- Thus for a typical BH-NS 
system with g ;^ 10 and R/M^s — 5, we have rcR > rm and general relativistic corrections to 
the orbital motion are expected to dominate over the Newtonian hydrodynamic effects discussed 
here. The situation is different for NS-NS binaries with q^l, where rcR ^ Tm typically. General 
relativistic corrections to our Newtonian treatment for NS-NS binaries are discussed in S5. 



2.5. Equilibrium Models for Two Unequal Masses 

General Darwin-Riemann models for two nonidentical stars can be similarly constructed 
with our energy variational method (see LRS2). In particular, models for two finite-size stars 
with different masses, radii, polytropic indices, adiabatic constants and spins can be constructed 
Such models can be used to describe binaries containing two nonidentical neutron stars. As in 
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Roche-Riemann binaries (§2.4), when the masses of the two stars are sufficiently difl'erent, a Roche 
Umit can exist prior to contact, providing a possibiHty of mass transfer. Here also, the Roche limit 
configuration is always secularly unstable, but can be dynamically stable or unstable depending on 
the mass ratio and compressibility. We return to these questions in §6, where wc model neutron 
star binaries as general Darwin-Riemann configurations in which the effective polytropic indices 
and radii of the two components are determined from a realistic equation of state and are functions 
of mass. 



3. ORBITAL EVOLUTION MODEL 

We now wish to study the orbital evolution of the binary models constructed in §2 in the 
presence of gravitational wave emission. As in §2, we consider binaries in circular orbits. This is 
probably justified for most systems at large separation since gravitational radiation itself tends to 
circularize an eccentric orbit. Near contact, however, if rclativistic effects are sufficiently strong, 
the eccentricity can actually grow again as the inspiral accelerates (Lincoln & Will 1988). This 
would be the case for two point masses approaching the last stable circular orbit allowed by 
general relativity. For neutron stars with finite radii, however, we expect hydrodynamic effects to 
become important before this relativity-induced eccentricity can grow significantly. We will return 
to this question in §5. Here, for simplicity, we assume that the orbit remains always circular, at 
least in some average sense. Wc determine the evolution of the average binary separation r as the 
system loses energy and angular momentum to gravitational radiation. Such an approach is clearly 
valid when the orbital decay time tr = \r/r\ is much larger than orbital period P = 2'k/Q.. But 
we adopt this approximation here even for the final phase of the orbital decay, when we find that 
the two timescales can become comparable. By doing so, we can study the transition from the 
secular orbital decay at large r to the dynamical coalescence at small r. We also assume that the 
internal structures of the two stars, and in particular, their shapes, assume the form for circular 
equilibrium configurations, i.e., we treat only r as a dynamical variable, assigning all internal 
degrees of freedom (pc, Ai, A2) to their equilibrium values. This approximation is valid as long 
as the internal dynamical time (the response time) tdyn ^ {R^/M)^/^ of the stars remains much 
smaller than the orbital decay time tr, a condition which is is well satisfied at large separation, 
and still marginally satisfied for two neutron stars near contact. 



3.1. Gravitational Radiation 

We calculate the emission of gravitational waves in the weak-field, slow-motion limit (see, 
e.g., Misner, Thorne &; Wheeler 1970). In this approximation, the rate of energy loss is given by 
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the usual quadrupole formula, (we set G = c = 1) 



where -iij is the reduced quadrupole tensor of the system 

-lij = J P (^^i^j - ■^^^'^y) (^^) 

For a binary system containing two stars of mass M and M' orbiting in the xy-plane, the only 
time-dependent components of the quadrupole tensor are 

-i^^ = [nr^ + (111 + I'n - I22 - I22)] ^ cos $ + constant, (42) 

-ixy = -lyx = [pr"^ + {hi + I'll - I22 - I22)] \ sin $ + constant, (43) 

-lyy = - [^ir^ + (111 + I'll - I22 - I22)] ^ sin ^> + constant, (44) 

where = MM' / (M + M') is the reduced mass, Iij and I-^ are the quadrupole moments of each 
star (cf. eq. Q), and we have defined an angle 

^> = 2 j VLdt + constant. (45) 

Expressions (^) and (^Tj) then give 



"'^ / GW 



1 H 2 (-^11 + ^'11 ~ -^22 - I22] 



For a circular orbit, the rate of angular momentum loss is given by (LRSl, Appendix D) 

-) =-{-) , (47) 
) GW ^ \dt J 

while the fluid circulation C is strictly conserved (Miller 1974). In expression (^), the second term 
in the bracket represents the correction to the point-mass result due to tidal effects. For large r, 
this is a small correction, of the order of Kn{R/r)^ (Clark 1977), but it can become as much as 
~ 40% near contact. 

In the quadrupole approximation, the wave amplitude hfj^ in the transverse-traceless (TT) 
gauge is given by 

hf;' = i'iIi\t-D), (48) 



where -ijj^ is the transverse projection of the reduced quadrupole moment, D is the distance 
between source and observer, t — Dis the retarded time, and a dot indicates a time derivative. For 
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wave propagation in the direction {9, </>) in spherical coordinates with orthonormal basis vectors 
e^, e^, and e^, the two basis polarization tensors are 

e+ = eg (g) eg - e^, (49) 
ex = eg (g) e^ + e^ (g) eg. (50) 

In this basis, hfj^ can be written 

1 ■• ■■ 2 

h^^ = /i+e+ + /ixex = ;^(-iee-%)e+ + ;^-ie0ex, (51) 

where the spherical components of the reduced quadrupole tensor are given by (we set (/> = 
without loss of generality) 

-igg = -ixxcos'^0, -I^^ = -iyy, -^^ = -ixyCos9. (52) 

Using equations (E2|) — (B5), (|5l|) and (|5^), we find the following expressions for the waveforms 



h+ = --r?2[^r2 + (/n -/22 -/22)](l + cos2 0)cos$, (53) 
hy, = --ir?2[^r2 + (In + 111 -122-/22)] cose sin (54) 

Note that the above derivation neglects the contribution from the orbital decay itself to 
the gravitational radiation. The correction to the wave amplitude due to the orbital decay is of 
order |r/(rr2)|, and the correction to the energy loss rate is of order r^/(rr2)^ ~ (M/r)^. This is 
smaller by a factor ~ (M/Rq)^ than the correction due to tidal effects that we have included in 
expression (p6|). 



3.2. Orbital Evolution at Large Separation 



For sufficiently large orbital separation, we expect the orbital decay to proceed quasi-statically 
(i.e., along an equilibrium sequence) with the rate of change r of the orbital separation given by 



dt ) GW \ dr 



(55) 



This expression must obviously break down when E(,q{r) is near a minimum, since it would 
otherwise predict that r ^ oo there. In this section we explore the effects of tidal interactions 
in the limit of large r, when expression (|55| ) applies. In particular, we calculate analytically the 
deviations from the point-mass behavior due to the finite-size effects. In §3.3, we will develop 
and implement a numerical formalism (ODE's) to calculate the orbital decay at smaller r, when 
equation ( ^5| ) no longer applies. 
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3.2.1. Point-Mass Results 



For binaries containing two point masses M and M' we have E^q = —MM'/{2r) and 
lij = I[j = in equation (^6|), so that equation (|5^) yields the famihar result 

— (50) 

5 r'^ 

where Mt = M + M' . The orbital evolution is obtained by integration, 

[r{tt = —^^Mf{T-t), (57) 

where T is the time at the end of the coalescence (r = 0). The frequency and phase of the 
gravitational waves are then easily obtained as 



fGw{t) = ^ = - 
vr vr 



256 ^Mt^^{T-t) 



3/8 

(58) 



/ T -t \ 

where is a constant and we have used = Mt/r"^, the Keplerian value. 

Of great importance for the detection of gravitational waves by laser interferometers is the 
number of orbits Norb (or the number of cycles of gravitational waves Now = "^^orb) in a given 
interval of wave frequency or binary separation (Cutler et al. 1993). This is obtained by integrating 

For point masses, this reduces to 

^^iS = T2^'^' = ^ 2/3 / f ^ .^,^ d{\nfGw). (61) 

1287r^Mf/^ 1927r^Mt^/^ {T^fGwr'^ 

Integrating equation (|6l| ) we find the number of orbits between and rj < rj, 

(0) _ 1 {rT-^T^ , . 

64vr MM'Ml^' ' ^ ^ 

Any deviation from this result leading to a change dN^rb ^ 0.06 over the detection interval 
(corresponding to a phase change 6^ = 2tt6Ngw = ^T^^^orb ^ ^/4) must be incorporated into the 
theoretical waveform templates used for signal extraction. 



We now discuss the deviations from the result ( |62[) caused by finite-size effects. We consider 
the three types of binary configurations introduced in §2: irrotational (C = 0) configurations 
(§3.2.2), configurations with C = constant/ (§3.2.3), and synchronized configurations (§3.2.3). 
For simplicity, we assume that M' is a point mass (i.e., we consider only Roche and Roche-Riemann 
binaries), but the generalization to two stars of finite size is straightforward in this case: one would 
simply add another contribution obtained by interchanging M and M' in the results given below. 
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3.2.2. Irrotational Configurations 

For nonspinning neutron stars (C = 0), the correction to iV^I'^ at large separation is entirely 
due to the tidal interaction. Since the tidally induced ellipticity of the star is ~ {Ro/rf, the tidal 
corrections to ft, {d£ /dt)GW , and dEeq/dr are all of order (Ro/r)^. Thus from equation (|60|), 
we expect a change 6{dNj^ll) oc {Ro/r)^dN^^^ oc r~'^/'^dr. A detailed calculation based on series 
expansions of the binary equilibrium equations at large r (see LRS2) gives for the total equilibrium 
energy 

„ , , MM' 3 

Eeq{r) = — + + ... + constant, 63 

(recall that g„ = Kn(l — n/5) and k„ is given by eq. |^]) and for the orbital angular frequency 

M + M' ( 9 M' Rl \ , 
''" = ^^[^ + r-'-M^ + --)- (64) 

The correction to {d£ / dt)GW due to the quadrupole moment of M (cf. eq. |^]) is 

//r^ 1 M r° 

Equation (|6^) then gives the change in N^j-iy for irrotational configurations^ 



We see clearly that this deviation from the point-mass result does not accumulate at large r. The 
term proportional to M~^M^ in expression (^) results from the change in Q. and E^q^ while 
the term containing the extra factor of Mt/M' comes from the increase in {d£ / dt)GW caused by 
the quadrupole moment of the star. 



3.2.3. Configurations with constant C 7^ 



For neutron stars with nonzero spin, the dominant effect at large r is the change in 
caused by the spin- induced quadrupole moment (Bildsten & Cutler 1992). Indeed, for 
sufficiently large r, tidal effects can be ignored, and the star can be modeled as an axisymmetric 
compressible Maclaurin spheroid, with In = I22 > ^33- The quadrupole interaction energy is 



^The numerical coefficients in our result do not agree with those given by Kochanek 1992 (cf. his eq. [5.4]). 
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-M'(/ii - /33)/(2r3) and we find (LRS2) that the total equihbrium energy can be written simply 



as 



MM' M' 
Eeq{r) = + ^(-^11 - -^33) + constant, 



while the orbital angular frequency is given by 

, M + M' 



0^ 



1 + 



2Mr2 



{hi - I33) 



(67) 



(68) 



The correction factor for {d£ /dt)GW is then simply [1 + 3(/ii — /33)/(2Mr^)]^. Therefore from 
equation (^0|) we have 



dNorb = dN. 



(0) 

orb 



1 



r(/: 



11 - -'33; 



2Mr2 

At large r this gives for the change in Ngj-b due to spin: 



1 + 



2Mr- 



11 — -^33 



-5/2 



SN 



105 (111 - 133) 



.1/2 1/2. 



orb 



2567r M^M' 



1/2 



(69) 



(70) 



This result agrees with equation (10) of Bildsten & Cutler (1992). Expressions for the quadrupole 
moments In and 133 of compressible Maclaurin spheroids can be found in LRSl. 

If the spin-induced eccentricity = 1 — al/af ^ 1, then the 0,s{e) relation for a rotating 
spheroid (e.g., LRSl, eq. [3.21]) can be expanded to give ~ (5g„/2)A2, with Cl^ = Vtl/{M/Rl). 
Thus we have: 



-^11 — -^33 _ of — al Kn 2 



MRl 

Therefore equation ( [70| ) becomes 

105 

K. 

5127r 



m 



2 



6N. 



(S) 



orb 



^i_j^2(^ 



1/2 l/2^ 



■MM' 



M^ 



1/2 



« 1). 



Since 5N^fl oc r^/^, we see that this deviation from the point mass result, in contrast to 5N^ 
does accumulate at large r. 



(71) 



(72) 



{I) 

orb ' 



3.2.4- Synchronized Configurations 

For the corotating case, the dominant effect is simply the added kinetic energy and angular 
momentum of the synchronized spin (Bildsten & Cutler 1992; Kochanek 1992). Expanding 
equation ( ^2| ) at large r, we obtain 

^ MM' 1 ,,„o/M + M'\ 

Eeq{r) = — + -KnMRl . + . . . + constaut. 73) 

2r b \ / 
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The corrections to {d£ /dt)GW ai^d ^ are of higher order in Rq/t, and we find from equation (|& 
that for synchronized spin, 

^^^orft 32^'^" ^^/2 

This result agrees with equation (3) of Bildsten & Cutler (1992) in the incompressible limit (where 
Kn = 1). Here also, since ^N^^^ oc r^/^, we see that the phase error accumulates at large r. 



3.3. Approach to Dynamical Instability 

Whenever E(.q{r) has a minimum at soinG T — ^iri! tliG orbital decay Cciniiot remain quasi-static 
as T^fi IS approached, since equation (^5|) would predict that T — > oo as T — > TjYi • This should not 
be too surprising since, as discussed in §2, the binary orbit becomes dynamically unstable for 
r < Tm- As the stability limit is approached, the radial infall velocity r can become much larger 
than predicted for two point masses (eq. For r < rm, the coalescence would proceed on a 

dynamical timescale even in the absence of energy and angular momentum losses. 

Let us first estimate the separation Vc > where equation (^) starts to break down. We 
can expand E(.q{r) around the minimum at r = as 

Eeq{r) = Eeq{rm) + —{r - r^f + . . . , e-M^/r-m- (75) 

Equation ( [5^ ) becomes invalid when the rate of increase of infall kinetic energy, becomes 
comparable to the rate of change of the equilibrium energy, i.e., is determined by 

/irr ~ ( — r, for r ~ Vc- (76) 
\ dr J 



Using equations ( ^5|) and (|75|), this reduces to 



dr 

where we have defined 

<5c = - rm)/rm < 1- (78) 
Using equations ( ^6| ) and (|77|), we then obtain 

Even for neutron stars, we have — 3i2, so that r^/M ~ 5 and is not very far from r^, 
typically 10% further out. 
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To properly calculate the orbital evolution for r < Tc, when the kinetic energy of the radial 
infall becomes important, we write the total energy of the system (not necessarily in equilibrium), 
as ^ 

£ = -fir^ + E{r,c/,;M,J,C), (80) 

where the second term is given by equation (^), and (a^ denotes the three variables {pc, Xi, A2). 
In writing down equation (pO|), we have implicitly ignored other possible contributions to the 
kinetic energy such as terms like Maf which are related to the change of the structure of the 
stars as the binary separation decreases. This is justified since the adjustment of the stellar shape 
takes place on the internal dynamical timescale tdyn ~ {R^/M)^^"^ of the star, while the orbital 
evolution of the binary usually takes place over a longer time tr- We also assume that the orbit 
remains quasi-circular, so that equations (^) and (^) hold. In doing so we neglect terms of order 
e^, where e is the eccentricity. Since ~ (r/fir)^, this is valid so long as tr ^ P. We find in §4 
that the combined condition tdyn <C P is well satisfied at large r, and remains marginally 
true even when r is close to r^. 



Taking the time derivative of equation (|8^), and recalling that gravitational radiation 
conserves C, we get 

o TP a p 

£ = f^rr + nj+—r + Y^—al (81) 

i * 

where we have used ft = dE/dJ. We now impose the assumption that the three parameters 
Pc5Ai,A2 specifying the internal structure of the stars take their equilibrium values. Accordingly, 
we have 

_ = 0, for a^ = (a:)e,. (82) 



Now since £ = [£)gw and J = {J)gw-: using equations ( |47| ) and (82), the evolution equation ( |8l| ) 
reduces to 

pr+—=0. (83) 



We now substitute expression (|20| ) for E and find on differentiating 

r-n\ + Ql^r = 0, with = — + -{2In - I22 - h3)- (84) 

In the last expression, the lij assume their equilibrium values as a function of r. We now express 



Q in terms of J using equations (|l5|)-(|lq). After some algebra we get 



n = j^^+F{r), (85) 



where 

Pir) = ^ It = ^lr' + ^MKn '^' . (86) 



^ _ 2 2 (Qi - 02; 

2 , 2' -fi — M' J. -i" f^n 2 , 2 

h af + 02 5 af + 02 



Equations ( |8^ ) — (|8^), together with equations (||) and ( |47| ) for j = {dJ/dt)GW can be 
integrated numerically given initial conditions at any separation rj such that (r^ — rm)/f'm ^ 
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We cast these equations into a set of first order ODE's for r, = r and J. We calculate initial 
values for r and r at r = rj from equation ( ]55| ) and then substitute r into equation ( p^ ) to obtain 
the initial value of J. 

Alternatively, we can eliminate J completely to obtain a single, third-order evolution equation 
for r. Taking the time derivative of equation ( |8^ and substituting expression (^) for J, we get: 

(fr rr ( 2r dIA f dF F dlt\ ,^ 2r dEeq 2r f d£\ 

l--:^)+2-M — + -T^)(^e,-f^) + --^r = -(-)^^, (87) 



dt^ r \ It dr J \dr It dr J It dr It \dt J 



where we have used QJ = E and 0,eqidJeq/dr) = dEeq/dr. A similar type of equation was derived 
by Lattimer & Schramm (1976) in their study of tidal disruption by black holes. For large r, the 
first three terms in equation ( |87| ) can be ignored since the acceleration of radial infall is small 
(note that {Qeg — ^) oc r from eq. |Q), so that equation ( |87[ ) reduces simply to r{dEeq/dr) = E, 
i.e., equation (|5| 



The derivation presented above assumed constant fluid circulation and is valid in the limit of 
zero viscosity. In the opposite limit, when viscosity is always acting on a timescale shorter than the 
orbital decay timescale, the binary remains synchronized throughout the evolution. The energy 
functional for Darwin (uniformly rotating) configurations should then be used in equation (^0|). 
We find that the orbital evolution in this case is still given by equations (|8^)-(|85|), but with 

F = 0, It = fir^ + 21 = nr^ + lMKn{al + al). (88) 

5 

We reemphasize that the assumption of uniform rotation must break down for r < r^ since 
viscosity will then drive the system toward lower energy, hence, away from, rather than towards, a 
synchronized state. Nevertheless, for simplicity, we calculate the orbital evolution in this limiting 
case based on the energy of Darwin configurations, even for r < r^- 



4. APPLICATIONS TO BINARY NEUTRON STARS 



We now apply the orbital evolution model developed in §3 to calculate the coalescence of 
two neutron stars. In §4.1, we first discuss how polytropes can best be used to approximate 
the internal structure of a neutron star. In §4.2, we examine the effects of the spin and tidal 
interaction on the gravitational waveforms at large separation, using the analytic results of §3.2. 
We then proceed to solve the orbital evolution equations derived in §3.3 and study the approach 
to dynamical instability at small separation (§4.3). 



4.1. Polytropic Models for Neutron Stars 
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The main parameter that enters the evolution equations of §3 is the ratio Rq/M^ which for 
neutron stars is determined from the nuclear equation of state (EOS). For the canonical neutron 
star mass M = 1.4M0, all EOS tabulated in Arnett & Bowers (1977) give Rq/M in the range of 
4-8. Small values of Ro/M correspond to a soft EOS, while large values correspond to a stiff EOS. 
Very soft EOS, such as that of Friedman &; Pandaripande (1981), appear have been ruled out by 
the constraint on the lower mass limit, 1.55M0, of the X-ray binary 4U0900-40 (Joss & Rappaport 
1984) as well as by the constraint derived from pulsar timing after glitches (Link, Epstein & Van 
Riper 1992). The most recent microscopic EOS constrained by nucleon scattering data and the 
binding of light nuclei are those of Wiringa, Fiks &: Fabrocini (1988, hereafter WFF; see also 
Baym 1991 for a review). For M = IAMq, the WFF EOS give Ro/M ~ 4.7 - 5.3, corresponding 
to Rq ~ 10.4-11.2 km. In addition, the radius is almost independent of the mass for M in the 
range of O.8M0 to I.SMq. Thus the value of Ro/M can be somewhat larger for smaller M. In this 
section we consider Ro/M = 5 and Ro/M = 8 as representative values. 

A polytrope is only an approximate parametrization for a real EOS. To find the approximate 
polytropic index n that mimics the structure of a real neutron star, we proceed as follows. For 
given M and Ro, we determine the ratio ///«, where / is the moment of inertia tabulated for a 
real EOS, and /„ is that of a uniform sphere with same M and Ro- Note that when we include 
general relativistic corrections, lu can be larger than 2MR'^/5. We calculate lu using equations 
(3.8)-(3.11) of Arnett & Bowers (1977). The resulting ratio is set equal to k„, from which the 
corresponding value of n is obtained. In Table 3, we list the results for different masses based on 
the EOS AV14+UVII of WFF. For the other two EOS given in WFF, the results are very similar. 
Typically, for M ~ 1.4M0, we find n ~ 0.5, highly incompressible. As M decreases, n increases 
and the configurations become more compressible. In the orbital evolution calculations presented 
below, we consider the representative values n = 0.5 and n = 1. We also give some results with 
n = 1.5 for comparison. 

When translating from dimensionless quantities to physical quantities such as a gravitational 
wave frequency in Hz, the values of both M and Ro/M must be specified. The masses of neutron 
stars have been determined for a number of binary radio pulsars as well as binary X-ray pulsars 
(see Thorsett et al. 1993 and references therein). All the measurements are consistent with a mass 
M = (1.35 =b O.1)M0. Hence we shall focus on this canonical value of M = I.4M0 when we quote 
actual wave frequencies. 



The initial neutron star spins also enter the calculation in the zero- viscosity case (cf. eq. [p6| ). 
In the four neutron star binaries known in our Galaxy, the radio pulsars have spin periods above 
30 ms, but much shorter pulsar periods ~ 1.5 ms have been observed in other systems. The 



minimum spin period corresponding to equation (37) is about 0.85 ms for M = 1.4 M0. Here we 
consider the representative values = 0, 0.1, 0.2, and 0.4, corresponding to no spin or spin 
periods of 4.8 ms, 2.4 ms, and 1.2 ms, respectively. We assume for simplicity that both stars have 
the same spin. 
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4.2. Orbital Evolution of Binary Neutron Stars at Large Separation 



The expressions derived in §3.2 can be evaluated easily for neutron star binaries. For 
definiteness, we consider the typical case where M = M' = 1.4 Mq, Rq/M = 5 and n = 0.5 (the 
corresponding value of Kn = 0.8148). We focus on the constant-C evolutions and assume that both 
stars have the same spin. The gravitational wave frequency few in I'eal units is given at large r by 

/ P \ -3/2 / X -3/2 

few = 5837 M{j — Hz, 



where Mi. 4 = M/(1.4Mq). For these parameters, the low-frequency band of interest to LIGO 
corresponds approximately to binary separations between = 70i?o {fGW,i = 10 Hz) and rj = 5Ro 
ifcwj = 522 Hz). When r ^ 5Ro, the analytic results of §3.2 become inaccurate. The orbital 
evolution for r < 5i?o is calculated in §4.3 using the method introduced in §3.3. The total number 
of cycles of gravitational radiation emitted by two point masses between and r/ is N§1. = 16098 
(eq. |6l). 

Using equation (^6[) we can calculate the accumulated change in Nqw = '^^^orb due to tidal 
effects in irrotational configurations. 



Ro\'/' fRo\'/^ 



-8.74 X 10-=Mf// (^) " - /Sg-.) , (90) 

where few is the wave frequency in Hz. This result is illustrated in Figure 3. Note that we 
have multiplied equation ( |66|) by 2 to account for the presence of two identical stars. The final 
change is 5NQ^{rf) ~ 0.30. Note that this change accumulates mainly at large few (small r). 
For few < 300 Hz, corresponding to r > 7.23-Roi or the first 16065 cycles, we find SNq^ <^ 0.1, 
while between few = 300 Hz and few = 522 Hz (the remaining 33 cycles) we find 5N^y^ ~ 0.2. 
Although these changes are small and can probably be neglected in the construction of theoretical 
low-frequency wave templates, they may nevertheless be detectable by advanced LIGO detectors. 

Now turn to the case where the stars are spinning. We assume that 17^ ^ 1 and use 
expression ([72|) (multiplied by 2 for two stars with the same spin) to get 



5/2 

5N}^> ~ -6.16 ( — ) 



\ 1/2 / ^ \ 1/2 



RoJ \Ro 



The total change is 5N^-^{rf) ~ 37.8 ~ 8.85/P^,, where 

Pms is the spin period in milliseconds. 

Clearly, this spin- induced change accumulates mainly at large r (see Fig. 3), and this effect can be 

( s) 

very important for rapidly spinning neutron stars: to get SNQ^^{rf) ^ 0.1, we need a spin period 
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Pms ^ 9. In agreement with Bildsten &; Cutler (1992), we conclude that finite-size effects for 
these rapidly spinning neutron stars are potentially very important in modeling the gravitational 
radiation waveforms, even at low frequency (large r). 

For synchronized binaries, the change 6Nq^^ calculated from expression (^) is much larger 
(by about two orders of magnitude) than either 5Nq^ or SNqy^. However, synchronized or 
nearly-synchronized configurations are particularly unlikely for binary neutron stars at large r. 
This is because the ratio of synchronization to orbital decay timescales increases with r (Kochanek 
1992). If viscosity plays any role at all during the coalescence, it is much more likely to be during 
the final phase. 



4.3. Orbital Evolution of Binary Neutron Stars Near Contact 



We calculate the orbital evolution at small r by integrating the evolution equations (|84|)~(8(:) 
numerically using a fourth-order Runge-Kutta method with adaptive stepsize to ensure accuracy. 
The integration is terminated at the separation r = rj where a contact configuration is reached 
along the equilibrium sequence. For two identical, slowly spinning neutron stars with mass 
M 0.7 and effective polytropic index n ^ 1 (cf. Table 3), we always have rj < r^, where is 
the separation where Eeq{r) is minimum (the stability limit, cf. §2). Since most of our assumptions 
in both §2 and §3 are only marginally valid when r < r^, the results for r/ < r < should be 
considered very approximate. 

Figure 4 shows the evolution of a system with n = 0.5 and Ro/M = 5. Constant-C evolutions 



with different values of the initial spin are shown. The point-mass results (eqs. [gg], [57|, and [plp 
are also shown for comparison, as well as the results for a corotating system. We see that, as the 
stability limit at r = is approached, the orbital decay rapidly accelerates and departs from 
the point-mass result, with r reaching typically about 10% of the orbital velocity near contact. 
As a result, it takes only one more orbit to complete the final decay from to rj. We see also 
in Figure 4 that the results for C ^ are not very different from those for C = 0. This comes 
about from two opposite effects: the initial spin tends to reduce the instability at small r (see §2, 
Table 2), while it can accelerate the orbital decay at larger r due to the spin- induced quadrupole 
interaction (see §3.2.3). By contrast, we see that the corotating evolution is very different from 
the constant-C evolution. 

Figure 5 compares the solutions for different n, all with Ro/M = 5. Clearly, for larger n 
(higher compressibility), the effects of the instability tend to be smaller. This is because a more 
compressible configuration is more centrally concentrated, and tidal effects only become important 
at smaller r. As a result both and — ?"/ (the "acceleration distance") are smaller (cf. 
Table 1). 
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During the evolution, the total angular momentum and energy of the system always decrease. 
This is in contrast to the equilibrium energy and angular momentum, which both increase for 
r < r^- This point is illustrated in Figure 6, where we show both the total energy of the system 
during the evolution £{r) and the equilibrium energy Eeq{r). At large r, ~ E^q, but when 

r ^ r-m, £ < Eeq. 

The rapid coalescence of the binary for r < rm can affect very significantly the gravitational 
radiation waveforms. In Figure 7, we show the wave amplitude /i+ seen by an observer along 
the rotation axis {6 = 0). In this case, hx differs from hj^ only by a constant phase of 7r/2 (cf. 



eqs. [53|-[£j]). Here we show the irrotational (C = 0) and corotating evolutions, and compare them 
with the point-mass result, for both n = 0.5 and n = 1, all with Ro/M = 5. We have set t = at 
r = 5Ro, and the various solutions all have the same phase at that point. We see that the phase 
and amplitude depart significantly from the point mass result, even when C = 0. The departures 
are largest for the corotating case (but recall that the phase error accumulates mainly at large r in 
that case; cf. §3.2.4). Figure 8 shows the number of cycles Nqw = "^^orb of the gravitational wave 
as a function of wave frequency during the evolution. The wave frequency in real units is given by 

= - = 4130 (^) ^^j^j^ Hz. (92) 

where VL is given by equation (|8^) (and remains only approximately equal to the equilibrium 
value, eq. |2^). We also see from Figures 6 and 7 that, for a given Rq/M, finite-size effects on the 
waveforms are more important for smaller n (stiffer EOS). 

In Table 4, we summarize our results for different cases. The first two blocks in the table 
correspond to constant-C solutions with = and Clg = 0.4. The last block corresponds 
to corotating evolutions. For every binary model, we list values of the stability limit r^, and 
the separation at contact, rj. In each case, we give the key orbital evolution parameters when 
Ro/M = 5 and Ro/M = 8: the radial velocity at r = r^, the radial velocity at r = rj, the 
maximum gravitational wave frequency f(^^^ = /gvk(^ ~ ^z); maximum wave amplitude 
hmax, and the number of cycles Nqw of gravitational radiation from r = 5Ro to r = rj. 



5. GENERAL RELATIVISTIC EFFECTS 



Our treatment so far has ignored Post-Newtonian (PN) effects other than the lowest-order 
dissipative effect corresponding to the emission of gravitational radiation according to the 
quadrupole formula (|40|). However, for the typical value of Ro/M = 5, other PN effects are likely 
to be important and can alter the orbits considerably. 

In the case of two point masses, these PN effects can make a circular orbit become unstable 
when the separation is smaller than some critical value ("inner-most stable orbit") ran- Kidder, 
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Will & Wiseman (1992) have recently obtained 

roB^ ~ 6Mt + 4iJ, (93) 

where Mt = M + M' is the total mas^. Their method includes PN effects up to order (w/c)^ 
in the treatment of the two-body interaction, but also incorporates test-particle effects in the 
Schwarzschild geometry exactly. Indeed, for a test particle orbiting a spherical black hole, 
expression (p^ ) reduces to the familiar Schwarzschild result rc/j = 6Mt. For two point masses 
with M = M', expression (^) gives the result[[| rc/j ~ 14M. Since — 3i?o typically, we see that 
rcR and are comparable for Ro/M = 5. 

For a test particle orbiting a spherical black hole of mass Mt, the total energy describing the 
orbital motion is given by 

'E\'^ f 2MA / J2 ^ 



where r is the Schwarzschild radial coordinate. The corresponding equilibrium energy for 
stationary circular orbits is 

'Ee,\' _ (r-2Mtf 
fi J r(r-3Mi)' ^ ^ 

which has a minimum at r = rcR = QMt- Although equation (|9^) is very different from our purely 
Newtonian expression (^), the two have one thing in common, namely the existence of a minimum 
marking the onset of instability. Thus we see that Newtonian tidal effects and relativistic effects 
both lead to the same qualitative result: the existence of a critical binary separation where a 
circular orbit becomes dynamically unstable. 

For two stars with comparable masses, a simple analytic expression such as (^) cannot be 
written since stationary circular orbits do not exist for a system radiating gravitational waves. 
Kidder, Will, & Wiseman avoided this problem by artificially turning off the radiation reaction 
terms in their PN equations of motion. In the same spirit, we adopt the following simple ansatz 
for the energy of the point-mass binary system 

^_2M,^ J^_2M. + 4./3 

The corresponding equilibrium energy for circular orbits is obtained by solving 
J M M' =0, which gives 

^ (r-2M,)^-(4/x/3)(3r/2-2M,) 

r{r-3Mt-2ij) ' ^ ' 




^Their original result was derived in harmonic coordinates. Here we ignore coordinate distinctions and fit their 
result to equation (^) as an approximation. 

^Note that this is very diflerent from the early estimate of Clark & Eardley (1977), who give ran ~ 6M. 
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and has a minimum at tgr given by equation (psl). In equation (p6|), the coefficient 4;u/3 was 
chosen to obtain this result. Note that when — > 0, equations ( |96[ ) and (p^) reduce to the exact 
results for a test particle, equations and (p5|). 



To estimate the combined effects of general relativity and the Newtonian tidal interactions for 
finite-size stars, we adopt the simple model introduced in Paper I. The equilibrium energy of the 
Newtonian fluid system discussed in §2 can be modeled approximately as 

eq 2r 2a r" ' ^ ' 

which has a minimum at r = r^- The parameters a and depend on the internal structure of 
the stars and the degree of synchronization. Their values are adjusted to obtain the best possible 
fit to the more accurate Ef^q calculated in §2. The total energy of the binary, not necessarily in 
equilibrium, is written as 

2\dt) 2ijr^ r 2a{a - l)r"' ^ ' 

For a stationary orbit (r = 0), the equilibrium condition {dE/dr)j^M,M' = yields the equilibrium 
energy given by equation (^). Taking the time derivative of equation (|99|), we obtain the orbital 
evolution equation 

J2 MM' MM'r^-^ , 
//r - — + +^ + = 0. (100) 

This equation, together with J = E/Q. = {fir'^ / J)E, can be solved for r{t). Equivalently, one can 
also eliminate J and derive an equation similar to (p 



(fr ^ 3/i fdr\ ( d'^r\ ^ 2 dEeq dr 2 fd£\ 



^dt^^ r [dt)[dt^j^r dr dt rKdtJcw' ^^^^^ 



The solutions of equation (101) reproduce all the essential features of the more accurate 
solutions obtained in §4. For example, for a typical case with = 2.8Ro, rj = 2.5Ro and a = 6, 
we get r(r„) = w,(r„) = -0.059(M/i?o) V2 and Vr{rf) = -0.12(M/i?o)^/^ for Ro/M = 5. These 
are close to the the typical values found in §4. 

We now attempt to incorporate the new relativistic effects discussed above into the model. We 
simply replace the Newtonian point-mass term —MM'/ (2r) in equation (|9^ ) by the corresponding 
GR term, expression (| 



MM' (rm 1 MM'r^'^ (rm , ^ 

Ee, = Ei^^ + ^ + eS""^ - M = ^ + eS^^ - (102) 

Of course, relativistic effects will also change the internal structures of the stars, changing the 
first term in the first equality of ( |102D , but this is a higher order effect. The new Eeq{r) now has 
a minimum at some which is larger than either the Newtonian or the purely relativistic 
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tgr- For example, for = 2.8Ro and Ro/M = 5, we find r!^ = SARq. To calculate the orbital 
evolution, we use the new Eeq{r) in equation ( |101| ) in place of Ei^\ In the example considered 
above, we obtain Vr{r'^) = -0M8{M/ Ro)^^^ and Vr{rf) = -0.25(M/i?o)V2^ i.e., the terminal 
velocity at contact is about a factor of two larger than when only Newtonian effects are considered. 
For Ro/M = 8, we have r'^ = S.ORo, and the terminal velocities are Vrir'^) = -0.028{M/ RoY^"^ 
and Vr{rf) = -0.13(M/i?o) ^2. We see that the general relativistic effects can be important, 
especially when the value of Ro/M is small. 

Our treatment of general relativistic effects is admittedly very crude. The main point we 
wish to emphasize here is that the Newtonian hydrodynamic effects discussed in this paper are 
likely to be at least as important as the relativistic corrections to the orbital motion for the final 
coalescence of neutron star binaries. When the two effects are combined, the final coalescence is 
likely to be even faster, and may assume a significant "head-on" character. 



6. POSSIBILITY OF MASS TRANSFER 



Our discussion so far has assumed that there is no mass loss from the system or mass transfer 
between the two stars. For nearly-equal-mass binaries, equilibrium configurations exists all the 
way down to contact (and even beyond; see Hachisu 1986), and there is no "Roche limit" in the 
conventional sense. When the two masses are sufficiently different, however, a Roche limit can 
exist and mass transfer becomes a possibility (see §2.4, 2.5). In particular, Clark & Eardley (1977) 
have suggested that stable mass transfer from the less massive neutron star to the more massive 
one can occur at the Roche limit. This stable mass transfer phase may last hundreds of orbital 
revolutions before the lighter star is tidally disrupted. It is accompanied by a secular increase of 
the orbital separation. Thus if stable mass transfer indeed occurs, a characteristic "reversed chirp" 
would be observed in the gravitational wave signal (Jaranowski & Krolak 1992). The problem has 
been reexamined more recently by Kochanek (1992) and Bildsten & Cutler (1992), who both find 
that very large mass transfer rates and extreme mass ratios are required for stable mass transfer 
between two neutron stars, making it rather unlikely. 

Our results suggest an additional problem with the Clark-Eardley scenario. Quite independent 
from the stability of the mass transfer itself, the Clark-Eardley scenario requires the existence 
of a dynamically stable Roche limit configuration. However, our studies of Roche-Riemann and 
Darwin-Riemann equilibrium configurations (LRSl, LRS2) indicate that this is almost impossible 
for objects as incompressible as neutron stars. Consider again the equilibrium energy curve E(.q{r) 
for a binary. The Roche limit at r = r^j^, when it exists, corresponds to the point where r has a 
minimum possible value for circular equilibrium prior to contact. However, before such a limit can 
be reached, the binary separation must pass through a value > where Eeq{r) is minimum. 
For equilibrium sequences with constant C, this minimum coincides with the dynamical stability 
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limit, and all binaries with r < Vm are dynamically unstable (cf. §2.4). Therefore, if viscosity can 
be neglected in the neutron star binary evolution, no dynamically stable Roche limit can exist, 
and stable mass transfer can never occur. 

In the opposite limit, when viscosity is so efficient that corotation can be maintained 
throughout the evolution, the mass parameter range which permits the existence of a stable Roche 
limit is very small. This is illustrated in Figure 9, which shows the different types of terminal 
binary neutron star configurations as a function of M and M'. The neutron star model is based 
on the WFF equation of state AV14+UVII, with the effective polytropic index given in Table 3 
(see §4.1). The diagrams are constructed from our general Darwin-Riemann equilibrium models 
allowing for two nonidentical polytropes (see LRS2 for details and for applications to other types 
of binary systems). For irrotational (C = 0) systems (Fig. 9a), the existence of a Roche limit 
requires a mass ratio q ^ 0.9 or g ^ 1.1 for M ~ 1.4M0. However, as discussed above, this Roche 
limit is always dynamically unstable. For corotating systems (Fig. 9b), the existence of a Roche 
limit requires a mass ratio g ^ 0.8 or g ;^ 1.2 for M ~ 1.4M0. Stable configurations can exist at 
the Roche limit but only if one of the stars has a very small mass, M ^ O.4M0. 

We conclude that even when the masses of the two neutron stars are different, stable mass 
transfer is nearly impossible. The final phase of the orbital decay is always a rapid coalescence, 
and a "reverse chirp" in the gravitational wave signal is not expected. 
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Fig. 1. — Equilibrium curves of total energy, total angular momentum and orbital angular velocity 
as a function of binary separation along various sequences with n = 0.5. Here 0^ = (2M/r^)^/^ 
is the Keplerian angular velocity. Solid line: irrotational (C = 0) Darwin-Riemann sequence; 
dotted line: Darwin-Riemann sequence with C = —0.0652, corresponding to an initial spin 
Cls = Qs/{M/Riy/'^ = 0.1 for both stars; short dashed line: C = -0.1304 {Qs = 0.2); long 
dashed line: C = —0.2607 (Clg = 0.4); light dotted-dashed line: corotating (Darwin) sequence 
ifR = 0). 

Fig. 2. — Variation of the three principal axes ai, 02 and 03 along the same sequences shown in 
Fig. 1. 

Fig. 3. — The accumulated change in the number of cycles of gravitational wave due to finite-size 
effects at large separation. The change is shown as a function of wave frequency fow (starting 
with zero change at few = 10 Hz). The two stars are assumed to be identical, with M = 1.4Mq, 
Rq/M = 5, and n = 0.5. The dotted line shows the tidal effects in irrotational configurations, 
5Nq^ (eq. [90]); the dashed lines show the spin-induced change 5Nq^ (eq. [91]), for VLg = 0.1 and 
fig = 0.2; the solid lines show the combined tidal and spin- induced effects {Sn'^qIt + <^^gw)- 

Fig. 4. — Terminal evolution of selected binary models shown in Fig. 1. Here (a) shows the infall 
radial velocity Vr = r, (b) shows the time (contact is reached at f = 0) and (c) shows the number 
of orbits, starting from r = 5Ro- In all cases the two stars arc identical with n = 0.5, Ro/M = 5. 
The solid line is for = (C = 0), the short-dashed line for = 0.2, the long-dashed line for 
Clg = 0.4, and the dotted-dashed line for corotating evolution. All curves terminate when contact 
is reached. For comparison, the dotted lines show the results for two point masses. 

Fig. 5. — Terminal evolution of selected binary models with different n. In all cases the two stars 
are identical and have Ro/M = 5. Here (a) shows the infall radial velocity and (b) shows the 
number of orbits starting from r = bRg. The solid line is for C = 0, n = 0.5, the short-dashed line 
for C = 0, n = 1, the long-dashed line for corotation and n = 0.5, and the dotted-dashed line for 
corotation and n = 1. The dotted line is the point-mass result. 

Fig. 6. — Total energy of the binary system as r decreases with time for identical configurations 
with n = 0.5 and Ro/M = 5. The solid line is for = 0, the short-dashed line for 0,s = 0.2, 
and the long-dashed line for corotation. The dotted lines are the corresponding equilibrium energy 
curves (as shown in Fig. 1). 

Fig. 7. — Gravitational radiation waveform just prior to contact. Here (a) is for Ro/M = 
5, n = 1, and (b) for Ro/M = 5, n = 0.5. The solid lines are for irrotational Darwin-Riemann 
binaries (C = 0), the dashed lines for corotating (Darwin) binaries, and the dotted lines for two 
point masses; the long-dashed lines show the envelopes of the waveforms. D is the distance from 
the binary to the observer located along the binary axis (0 = 0). 
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Fig. 8. — Number of gravitational wave cycles Nqw = '^^orb as a function of the wave frequency 
few just prior to contact. Here Ro/M = 5 and M = IAMq. The solid line is for an irrotational 
binary with n = 0.5, the short-dashed line for an irrotational binary with n = 1, the long-dashed 
line for a corotating binary with n = 0.5, and the dotted-dashed line for a corotating binary with 
n = l. The light dotted line shows the point-mass result for comparison. 



Fig. 9. — Final fates of coalescing neutron star binaries. Each point in the diagram corresponds 
to a sequence of binary equilibrium configurations with different values of M and M' . The various 
regions indicate the possible terminal configurations. Neutron star models are based on the WFF 
equation of state. Irrotational systems are shown in (a), corotating systems in (b). In (a), no 
binary equilibrium sequence in the hatched regions has a Roche limit. Instead, a contact equilibrium 
configuration exists. In the wide-hatched region, this contact configuration is dynamically unstable, 
while it is dynamically stable in the narrow-hatched region. A binary sequence in the unshaded 
region does terminate at a Roche limit, but the Roche limit configuration is always dynamically 
unstable. In (b), no binary sequence in the shaded regions inside the solid lines has a Roche 
limit. The contact configuration is dynamically unstable in the wide-hatched region, while it is 
dynamically stable in the narrow-hatched region. Binary sequences outside the solid lines have 
Roche limits. A binary in the unshaded region encounters a dynamical instability before the 
Roche limit; this instability occurs after the Roche limit (along the second, unphysical branch of 
the equilibrium sequence, cf. §2.4) in the cross-hatched region. A binary in the blackened region 
encounters a Roche limit but no dynamical instability. Only in the cross-hatched and blackened 
regions can a dynamically stable binary exist at the Roche limit. Binaries in all regions in (b) 
encounter a secular instability (energy minimum) prior to any other critical point, except those in 
the small unshaded portion in the lower left corner. 
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Compressible Darwin-Riemann Sequences With C = 
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" Here Q., J and E are defined in eqs. (34); R and Ro are defined by eqs. (26)-(27); 
Ts is the "spin" kinetic energy (cf. eq. [14]) and W is given in eq. (3). 
* marks the dynamical stability limit. 



TABLE 2 

Compressible Darwin-Riemann Sequences With Constant C " 
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" Here Cl, J and E are defined in eq. (34); R and Ro are defined by eqs. (26)-(27); 
Ts is the "spin" kinetic energy (cf. eq. [14]) and W is given in eq. (3). 
* marks the dynamical stability limit. 



TABLE 3 

PoLYTROPic Model for Neutron Stars 



M{Mq) 


R (km) 


/ (10"^ g cm^) 
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" Based on the AV14+UVII EOS of WFF; I is the moment of inertia; 7„ is 

the moment of inertia of a uniform rclativistic sphere with the same mass 
M and radius R in Schwarzchild coordinates; fv„ = I/Iu- 



TABLE 4 
Binary Evolution Parameters " 



n 






Vm 




J GW \ / '''max 


Now 



fis = 0.5 2.76 2.62 



Ro/M = 5 








0.062 


0.12 


1302 


0.329 


15.1 


R„/M = 8 








0.040 


0.097 


653 


0.212 


47.6 




1.0 


2.49 


2.44 












R„/M = 5 








0.072 


0.090 


1518 


0.360 


16.8 


R„/M = 8 








0.044 


0.063 


761 


0.232 


53.5 




1.5 




2.29 












J?o/M = 5 










0.068 


1692 


0.375 


18.0 


Ro/M = 8 










0.033 


845 


0.240 


57.7 


fis = 0.4 


0.5 


2.67 


2.61 












R^/M = 5 








0.077 


0.099 


1387 


0.364 


14.6 


Ro/M = 8 








0.051 


0.072 


696 


0.235 


45.9 




1.0 




2.48 












J?o/M = 5 










0.075 


1515 


0.370 


16.1 


Ro/M = 8 










0.043 


758 


0.237 


51.2 


Corotation 


0.5 


2.99 


2.63 












R^/M = 5 








0.066 


0.17 


1261 


0.308 


10.1 


Jio/M = 8 








0.040 


0.14 


633 


0.198 


30.3 




1.0 


2.76 


2.52 












Ro/M = 5 








0.074 


0.14 


1391 


0.323 


12.0 


7?o/M = 8 








0.045 


0.11 


698 


0.208 


36.9 


" Here r = r/Ro, v 


= v/{M/Ro) 


is 


the point where Eeq{r 


) is minimum; 


r/ is the 



separation at contact; Vm is the value of the radial velocity at rm, and Vf is the velocity 

at r/; /q™^' = Qf /tv is the maximum frequency of the gravitational wave just prior to 
contact for M — 1.4M0, and hmax = hmaxM/D is the maximum wave amplitude for 
an observer at a distance D situated along the binary axis; Now is the number of cycles 
of gravitational wave from r = 5Jio to r = r/. 



